Solitons in the Yakushevich model of DNA beyond the contact approximation 
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The Yakushevich model of DNA torsion dynamics supports soliton solutions, which are supposed 
to be of special interest for DNA transcription. In the discussion of the model, one usually adopts 
the approximation to —* 0, where to is a parameter related to the equilibrium distance between 
bases in a Watson-Crick pair. Here we analyze the Yakushevich model without to — > 0. The model 
still supports soliton solutions indexed by two winding numbers (n, m); we discuss in detail the 
fundamental solitons, corresponding to winding numbers (1,0) and (0,1) respectively. 



Introduction 

Following the pioneering paper and proposal by Eng- 
lander, Kallenbach, Heeger, Krunihansl and Litwin 9], a 
number of authors considered simple idealized models for 
the roto/torsional dynamics of DNA one, with the aim of 
describing nonlinear solitonic excitations which - accord- 
ing to the ideas of Englander et al. - would be related to 
the transcription bubble which travels along with RNA- 
Polymerase in the DNA transcription process, and would 
thus play a functional role in this process. 

These models are based on modelling the DNA 
molecule as a double chain of coupled pendulums; the 
relevant nonlinear excitations would then be (topologi- 
cal and dynamical) solitons like those of the sine-Gordon 
equation. For a review of the approaches in this direction, 
see [23; for general properties of DNA and its functions, 
see e.g. [g- 

It should be noted that in a related approach - but 
in a different direction - the stretching motions of the 
DNA molecule (related to DNA denaturation) have also 
been studied by nonlinear models, see in particular 
the Peward-Bishop model [la| and extensions thereof 
[H Hi 113 ■ In this note we will focus on roto/torsional 
dynamics, and thus we will not deal with Peyrard-Bishop 
like models. 

Coming back to rotational dynamics, a very interesting 
and quite successful model (also called Y-model) was put 
forward by prof. L.V. Yakushevich 0. See [H 113,123 
for further results and extensions. 

In the Yakushevich model, the DNA molecule is con- 
sidered homogeneous (i.e. all nucleotides are considered 
as identical), and the state of each nucleotide is described 
by a single degree of freedom; in fact, a rotation angle 
of the base belonging to the nucleotide around the Ci 
atom in the sugar-phosphate backbone. Each nucleotide 
(or base) is represented as a disk. 

Interaction between successive nucleotides on each 
DNA helix is via a harmonic potential, while the intrapair 
interaction (i.e. interaction between bases in a Watson- 
Crick pair, and through these between the correspond- 
ing nucleotides) is modelled by a potential which albeit 
harmonic in the distance, becomes anharmonic when de- 
scribed in terms of the relevant degrees of freedom, i.e. 



rotation angles. More precisely, with I the distance be- 
tween relevant points Bi and B2 (see fig. 2) of the disks 
representing nucleotides in the Yakushevich model, and 
£0 the distance between the points Bi and B^ in the equi- 
librium (i.e. the B-DNA) configuration, the intrapair po- 
tential is Vo = (l/2)-K^p(^ — £0)^1 with Kp a dimensional 
constant. 

In the standard Yakushevich model, one considers the 
approximation (.^ = Q, which leads to a number of com- 
putational simplification; we call this the contact approx- 
imation^ as it corresponds to having contact between the 
disks representing nucleotides at same site on the two 
chains. However, as pointed out by Gonzalez and Martin- 
Landrove [I3I, ^0 = is a singular case in that the de- 
scription one thus obtains is not structurally stable: as 
soon as we consider (.q ^ 0, certain qualitative features 
of the model dynamics are changed. 

In this note we will analyze the Yakushevich model 
beyond the contact approximation, i.e. with a nonzero 
value for the parameter ^0; we focus in particular on trav- 
elling wave solutions and soliton-like excitations, as they 
are the most important objects in the approach of Eng- 
lander et al. 9]. 

It turns out that in this case soliton solutions are still 
present, albeit the simple form obtained with the contact 
approximation is replaced by an expression involving el- 
liptic integrals. The qualitative form of soliton solutions 
is little changed with respect to the standard Y-model; 
as for their width, it is very moderately increased (see 
fig. 5), still remaining in the same order of magnitude. 
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Dynamics (M'^ x D"^)". 



I. THE MODEL 

In the Yakushevich model the DNA double chain is 
considered to be infinite, and all bases are considered as 
having the same physical characteristics; it is thus an 





FIG. 1: Left; A base pair in the Yakushevich model. The 
only moving elements in nucleotides are the bases; these are 
represented as identical disks of radius r, which can rotate 
around their centers. The disks centers lie at a distance A 
from the double helix axis, and the rotation angles (positive 
in counterclockwise direction) are 6i, 62 respectively. The 
distances A and r are related via A = r + ^o/2, and la is 
the distance between the disks, Iq = 2{A — r). Right: disks 
represent bases and are located on the double helix; two disks 
as those marked in black interact via the helicoidal terms 
considered in the appendix. 



are the stacking [Us) term and the pairing one {Up). Thus 
U = Us + Up. (There is a third term Ut if we consider 
the so-called "helicoidal" version of the model j3, \W{ : 
introducing Uh leads to quahtative changes in the small 
amplitude regime. However, the additional term Uh is 
not relevant to fully nonlinear dynamics, and thus we 
will not consider it. See the appendix for small amplitude 
dynamics, where it matters to introduce this term.) 

The first term corresponds to a harmonic potentials 
(the choice of harmonic potentials for these term is com- 
mon to nearly all DNA models |15ll20l |) and depend on 
a dimensional coupling constant K^. It describes back- 
bone torsion and stacking interactions between bases at 
nearby sites on the same chain: in our notation it is given 
by 



Us 



2 



EEO 



i(a) 



Q(a) 



(1-2) 



As for the Up (pairing) term, it describes the nonlinear 
interactions between bases in a pair, i.e. 




FIG. 2: Further detail of the base pair modelling in the Yaku- 
shevich model, showing in particular the points Si and B2 
whose spatial distance is I. For graphics convenience, I and 
Ifj have been denoted as L and Lq in this plot. 



"ideal" DNA model according to Yakushevich's classifi- 
cation [1II23. 

Moreover, each nucleotide is considered as a single unit 
and its state described in terms of an angular variable. 
More precisely, each nucleotide is modelled as a rigid disk 
which can rotate around an axis and has a moment of in- 
ertia / around this axis, see fig.l. The rotation of the 
nucleotide at site i e Z on the chain a (a = 1,2) is de- 
scribed by an angle d^ ; we orient all angles in counter- 
clockwise direction. When we are referring to a specific 
base pair, we write simply dx, i?2 for -di , i?^ . 

The model is described by a Lagrangian L = T — U. 
The kinetic energy is 



/ 



r = EEi(C0 



q(a)A2 



(1.1) 



The potential energy U is the sum of two terms; these 



t/„ 



E 



Vb(i9 



(1) ^(2) 



(1.3) 



The potential Vq is not harmonic in terms of the angles 



^, 



W. 



here lies the nonlinearity of the dynamics and hence 



the hearth of the Yakushevich model. It is assumed that 



Vo = (1/2) Kp{£ - io)' 



(1.4) 



where i is the distance between the atoms which are 
bridged by the H bonds - these are represented in the 
model by reference points on the border of the disk repre- 
senting the bases, and more precisely by the point which 
lies nearer to the double helix axis for i? = - and io is 
this distance in the equilibrium configuration described 

by 'di' = 0. We also write r + £o/2 — A for the distance 
between the center of the disks representing nucleotides 
and the axis of the double helix. 

This harmonic potential becomes anharmonic in terms 
of the angles "di" . Indeed, with simple algebra (and writ- 
ing -da for "iJj ) wc have 



{£l + 4£or + 
+2r^ cos(i3i 



cos l?9 



6r2)-2r(4 + 2r)[cosi?i 

- ^^2) . 

(1.5) 

In the standard Yakushevich model, one considers the 

contact approximation Iq — 0, which entails A — r; with 

this, and omitting a constant term which plays no role. 



u 



Kpv' 



0(2) ^ 



g(l) 



q(2) 



> -^Y )-'2{co^^l +co^^i )]■ We 
for further detail on the standard 



y^.-[cos(i^ 
refer e.g. to [ij, [20^ 
Yakushevich model. 

In this note we will not adopt the contact approxima- 
tion £0 ~* Oi E^nd instead consider the potential energy 
(1.3) with Vq given by (1.4), (1.5). 

Let us provide an explicit expression for Vb(^i,i?2)- 
With the notation introduced above, i'^ = (2A— rcost^i — 



r cos ^2)'^ + {r sin -di+r sin i?2 ) ^ ■ Introducing the adimen- 
sional parameter A = r/A, this yields 



Vo= il/2)Kp{i-eo)^ 



2KpA^ 



(1-A)- Vp(^i,^2) 
where we have written for short 



pi^U^2) ■■= (1-A(C0S^1+C0S^2) + 

+ (AV2)[l + cos(79i-792)]) 



(1.6) 



(1.6') 



Needless to say, Vb('!^i, i?2) > for all values of ^1, ■(?2, the 
equality being satisfied only for the equilibrium position 
^1 = l?2 = 0. 

It should be stressed that this minimum of Vq is non- 
quadratic, as remarked by 13]; indeed, expanding around 
the equilibrium t9i = ??2 = 0, we get 

Fo(£i?i,ei?2) - {Kp/2){r[r{{},-^2?+ 

The equations of motion - i.e. the Euler-Lagrange 
equations arising from the modified Lagrangian - are 



0n^ = K. 



2# 



■d 






(1.8) 



It is convenient, as in the standard Yakushevich model, 
to pass to variables 



V'n 



^i'^ + ^i? 



Xn 



^(l)_^(2) 



(1.9) 



these correspond to ^„ = ■)/)„ + Xn and i?„ — 'i'n " Xn- 
In terms of the il),x variables the eqs. (1.8) read 

I^n = Ks (V-n+i - 2^^ + ^„_i) - (1/2) {dV/d^k) , 
IXn = Ks (x„+i - 2x„ + Xn-i) - (1/2) (dV/dx) . 

(1.10) 
Here we have denoted by F = V{^p, x) the expression of 
Vb(^i,i92) in terms of the new variables; that is. 



V{-ip,x) := Voiip + Xi'ip-X.) 



(1.11) 



II. PHYSICAL VALUES OF THE PARAMETERS 

Several dimensional parameters appear in our model 
Lagrangian and hence in the Euler-Lagrange equations 
of motion. We will discuss the model and its predictions 
for generic values, but specific values apply to the DNA 
molecule; for the geometric ones we will refer for dcfinite- 
ness to B-DNA. 

The parameter A represents the distance from the cen- 
ter of rotation of the disks representing nucleotides to 
the axis of the double helix. Taking this center at the Ci 
atoms (as in the derivation of parameters in the standard 
Yakushevich model), we get 2 A = ll.lAfor AT pairs and 



2A = lO.SAfor GC pairs ^. We will take the interme- 
diate value A = 5.48A. (Had we taken the center of ro- 
tation on the phosphodiester chain, A ~ lOAwould have 
resulted ^). 

The choice of the Ci atom as center of rotation for 
the nucleotides is not only conformal to discussion of 
the standard Yakushevich model [l^, lO . l20l |. but also 
physically sound: indeed the phosphodiester chain in the 
backbone is a very flexible polymer (as also confirmed 
by the success of Poland-Scheraga type models based on 
such flexibility 0, [lJ|)i while the complex formed by the 
sugar ring and the attached base can be seen in a flrst 
approximation as a rigid body |23 |. 

As for r and £0, these satisfy 2 A = 2r + £0; thus r is 
obtained as r = A — £q/2. Here £q is the length of the 
H bond bridging bases in a pair, while r is the radius 
of disks representing bases. The physical meaning of r 
would be the distance from the center of rotation for disks 
(the Ci atom in the sugar ring) to the atoms bridged by 
the intrapair H bonds. 

This distance is quite different from one base to the 
other, and also different for different H-bonded atoms in 
the same base. Thus we prefer to set the parameter of our 
idealized model in terms on the much more uniform value 
oi £0, which is £q ~ 2.9±0.lAfor the different H bonds in 
Watson-Crick pairs (see sect. 7. 2 of 17]). We will adopt 
the value £1^ = 2.9A. With these choices, r ~ 4.0A, and 
hence A = r/A = [1 - £o/{2A)] ^ 0.74 ~ 3/4. 

Finally, 5 represents the interbase distance along the 
double helix axis; in B-DNA this is (5 = 3.4A. 

After discussing the geometrical parameters, let us now 
come to the dynamical ones. First of all, we have the 
moment of inertia /; this is rather different for different 
bases (detailed values are given e.g. in [221 )• Taking an 
average over different values, we will adopt as in 11] the 
value / = 3 • lO^'^^cm^g. As for the coupling constants 
Ks and KpV^ (and Kh considered in the appendix), we 
will also adopt the values given in [IJ], i.e. Kg = 0.13eV, 
Kpr^ = 0.025eV, R\ = 0.009eV. 

In our discussion of soliton solutions and their width, 
we used the parameters /i and B; the first of these de- 
pends on the speed of the soliton, which is a free parame- 
ter (provide it is smaller than the limiting value y'Kj/IS, 
see sect. 3) in the Yakushevich model. For w = 0, we get 
fi = —Kg5^. The parameter B~^ = b6 is defined in terms 
of fj, by (5.8) below. With our choices for the dimensional 
parameters, and fi as above, we get b ~ (3/8)2.28 ~ 0.86, 
and hence B~^ ~ 0.865, B ~ 0MA~^. 



III. CONTINUUM APPROXIMATION 

In discussing (1.10), it is convenient to promote the 
arrays {^pnit)} and {xnit)} {n & Z) to fields i^{x,t) and 
x(x, t) (no confusion should be possible between old de- 
pendent variables and fields), the correspondence being 
given by 



i'nit) ^ Tp{nS,t) , Xnit) ^ xinS,t) 



(3.1) 



Here S is the spacing between successive sites of each 
chain; in B-DNA we have S = 3.4A. (It is of course also 
possible to pass to adimensional units in the spatial vari- 
able, i.e. set ^ = x/6, and consider tp{S,, t), x{£., t) so that 
e.g. Vn(i) = V'(C,*) for C = n.) 
With this, (1.10) reads 

Ii>tt{x, t) = Ks {ij{x + (5, t) - 2V'(x, t) + ij{x - 5, t)) + 

-{ll2){dV/dtl^){x,t), 
IXttix, t) = K, {x{x + S,t)- 2x(x, t) + x{x - <5, t)) + 

~{i/2){dv/dx){x,t). 

(3.2) 
If now we assume that ^(a;,i) and x(a;,t) vary slowly 
in space compared with the length scale set by lattice 
spacing, we can write 

il^ix ±S,t) = i/'(x, t) ± Sip^{x, t) + {5'^ /2)ip^^{x, t) , 
X{x ±6,t) = x{x, t) ± 6xAx, t) + {6^/2)x^,{x, t) . 

(3.3) 
Inserting these into (3.2), we obtain the equations of 
motion for the Yakushevich model in the continuum ap- 
proximation: 

I^tt{x,t)^ KsS^^P^,{x,t)-{l/2){dV/di;){x,t) , 
IXtt{x,t) - KsS^xUx^t) - (1/2) [dV/dx) ix,t) . 

(3.4) 
We omit from now on to specify at which point functions 
should be evaluated, as we got a local formulation. We 
will also consider, where appropriate, ^i and i?2 as fields. 
The PDEs (3.4) should be supplemented with a side 
condition specifying the function space to which accept- 
able solutions belong. The physically natural condition 
is that of finite energy; that is, we should require that 
the integral 



-HCXD 



1, 



- {li^f + x\) + KMl + xD) + 2^(V',X) 



Ax 

"(3.5) 

is finite; if this condition is satisfied at t = 0, it will be 
so for any t. 

It should be mentioned that the Yakushevich equa- 
tions (3.4) correspond, in the ^o = approximation, to 
classical ones when only one field is nonzero: indeed for 
X = they reduce to the sine-Gordon equation, while for 
-0 = one gets the so-called "double sine-Gordon" equa- 
tion, which appears in many physical contexts [J, Q . 



IV. TRAVELLING WAVE SOLUTIONS 

Next we focus on travelling wave solution for (3.4). 
That is, we restrict (3.4) to a space of functions 

"0(3;, t) = ip(x — vt) , x{x,t) = rj{x — vt) (4.1) 

(we will further restrict this in the following, in order to 
take into account the finite energy condition). We will 
also write simply z :— x~vt, and introduce the parameter 



Note that /i could be negative; this will actually be the 
interesting case. 

With (4.1) and (4.2), the (3.4) reduce to two DDEs, 
i.e. 



^" = ~[l/{2^^)]{^V{^,v)/^^) 

if = -[l/{2^^)]{^v{^,v)/^7J) 



(4.3) 



These describe the motion (in the "time" z) of a point 
particle of unit mass in the effective potential 

W{^,Tj) ■■= (2^)-iy(^,,7). (4.4) 

The conservation of energy reads then 

liiv'f + Wf) + W{^,rj) = E. (4.5) 

The finite energy condition (3.5) implies, in terms of 
<y9(z), 77(2;), that 

lim if (z) = lim ry (z) = ; (4.6') 

Z — ^±00 Z — ^±C30 

moreover, the functions (p and rj themselves should go to 
a point of minimum for the potential V. 

If /i > 0, the minima of V are the same as the min- 
ima of W; but if /i < 0, then minima of V are the same 
as the maxima of W. As it is impossible that nontriv- 
ial motions reach asymptotically in time (that is, in z) a 
minimum of the effective potential - while they can reach 
asymptotically a maximum if they have exactly the cor- 
rect energy - in order to have travelling wave solutions 
satisfying (4.6') and going to minima of V for z — *■ ±00 
we need 



M 



< 



(4.7) 



we assume this from now on. Note that /i < implies 
there is a maximum speed for travelling waves: 



\v\ < VKs/IS 



(4.8) 



The minima of V, i.e. the maxima of W, are for 
(1^1,^2) = (2gi7r,2927r); writing n = {qi + q2)/2, m = 
(91 ^ 92)72, these correspond to {(p,r]) = {2mr,2mTT). 
Thus the finite energy condition requires (with obvious 
notation) 



lim (p{z) = 2n±TT 

: — >±oo 



lim T]{z) ^ 2m±7T . (4.6") 



^ := Iv^ 



Ks6' 



(4.2) 



We can and will always take n_ = to_ — with no loss 
of generality; we will hence write simply n, m for n+, m^. 
These satisfy n — mG Z (in addition to 2n E Z, 2m G Z). 

In terms of the dynamical system describing the evo- 
lution in "time" z in the potential W, the solutions sat- 
isfying (4.6) represent heteroclinic solutions connecting 
the point (0, 0) at z = — cxd with the point (2?™, 2TTm) at 
z — +00. It is thus no surprise that the analytic deter- 
mination of such solutions is in general impossible. 

On the other hand, the solutions with indices (1, 0) and 
(0,1) can be determined explicitly, as shown in the next 
section. 



V. SPECIAL SOLUTIONS 

The solution with indices (1,0) and (0,1) are special 
in that they require that only one of the two fields varies. 
That is, the (1,0) solution will have ri{z) = 0; and the 
(0, 1) solution will have (p(z) = 0. Thus, they corre- 
spond to one-dimensional motions in the effective poten- 
tial W{(p, rj), and as such they can be exactly integrated. 

Note that these are immediately taken back to a de- 
scription in terms of the original angles "di: indeed, by 
(1.9), for the (1,0) solution we will have i?! = i?2 = <P, 
while for the (0,1) solution it results di = —'d2 = V- 



A. The (1,0) solution 

Setting ry = 0, the effective potential reduces to 

P(^) :=. W{^,0) ^ {2fi)-'V{^,0); (5.1) 

note that /i < implies P{ip) < for all ip. The first of 
(4.3) reduces to ip" = —[dP{(p)/d(p], and the conserva- 
tion of energy (4.5) reads 



P((P) 



(l/2)(^')' + Pi^) 



E . 



(5.2) 



By construction the solutions satisfying the side condi- 
tions (4.6) correspond to E = P(0), as 0(— oo) — 0; more- 
over, again by construction, 

P(0) = W{Q,Q) = {2n)-^V{Q,Q) = 0. (5.3) 

Thus, the separable equation (5.2) yields 



|^V-3?P 



(5.4) 



(In the (1,0) solution, ip' > 0; i.e. we have the positive 
determination of the root at all z.) 

The expression for P is readily obtained once we note 
that rj = {) means ^i = 1^2, see (1.9). With this, the 
expression (1.6) for Vb('!?i,^2) gets simplified. We write, 
for ease of notation, 



h{ip) := ^l - 2Acos(p + A2 , 
/3(^) :- ((A-1) + h{^)) ; 



(5.5) 



note that h{(p) > and /3{(p) > for all ip (equality 
applying only at ip — 2fc7r) for < A < 1. It results with 
standard algebra that 



Vo{d,d) = 2KpA^ (3^^). 



(5.6) 



It follows immediately from (5.6) and (4.4), (5.4) that ip 
is obtained by integrating 



^ B dz 



with B a constant, 



dip 

W) 



B:= AJ~2Kp/^i =. Aj2Kp/\^i 



(5.7) 



(5.8) 




f(q>) 



q> 



FIG. 3; Upper plot: the function P{ip), representing the 
"speed" Atp/Az (in terms of the effective dynamics descrip- 
tion provided by (5.4)) of the (1,0) soliton in adimensional 
units (see text and (5.7)), for A = 3/4. Lower plot: the func- 
tion f{ip) = J[l//3(<p)]dv?,to be inverted in order to get the 
(1,0) soliton solution (see 5.11)), again for A = 3/4; here co 
has been chosen according to (5.10), so that /(tt) = 0. 



That is, I3{ip) represents the rate of variation with z of 
the soliton solution in adimensional units, set by (5.8). 

The integral f{ip) of the left hand side of (5.7) is given 
explicitly by 



/(^) - Co - [2A]-i X [{l-X)E{^/2,<7) 
-[(l + A)V(l-A)]F(^/2,^) + 
+ [l-\ + h{ip})] cot((^/2)] . 



(5.9) 



Here F and E are the elliptic integrals of first and 
second kind respectively, defined as F(a;, cr) = J„ (1 — 

asm^ ey'^/'^de imd E{x,a) = ^^{l - asm^ ef^de. The 
complete elliptic integrals are JC{a) — F(7r/2, cr) and 
£{<j) = E(Tr/2,(T). In (5.9) we have moreover a = 
— 4A/(1 — A)^, and cq is the integration constant. The 
latter can and will be chosen as 



Co 



(1-A)2g(a) - (1 + A)2/C((t) 
2A(1-A) 



(5.10) 



with /C and £ the complete elliptic integrals of first and 
second kind respectively; in this way f{ip) is antisymmet- 
ric with respect to (ys = 0. The function f{ip) is singular 
at if — 2m:, as seen in fig. 2. 

Needless to say, integrating also the right hand side 
of (5.7), we get f{ip) = B{z - zq). With zq = (this 
integration constant can be absorbed in cq), this yields 
finally for the (1,0) solution (shown graphically in fig. 3) 



ip = f-\Bz) 



(5.11) 
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FIG. 4: The (1,0) soliton for the Yakushevich model without 
the contact approximation, see (5.11), with A = 3/4 (solid 
lines); and for the Yakushevich model with the contact ap- 
proximation (dotted lines). Here ip is in units of tt, and z in 
units of the distance 5 between successive base pairs. 



FIG. 5: The (0,1) soliton for the Yakushevich model with- 
out the contact approximation, see (5.5), with A = 3/4 (solid 
lines); and for the Yakushevich model with the contact ap- 
proximation (dotted lines); units as in fig. 3. 



B. The (0,1) solution 

For the (0,1) solution we can set (p = 0, which means 
'(?2 = — '^i- With this, and writing again r ~ XA, 



(£-4)' = 4^2 A^ (1-COS77)' 
The effective potential is hence given by 



(5.12) 



Q(r/)=iy(0,r;) 



ViO,rj) _ {AX)^Kp 



2^ 



M 



(1 — cosr;)^ 



(5.13) 

Note that again /i < implies Q{ri) < for aU rj, and 
that Q(0) — 0. Conservation of energy provides in this 
case dr]/dz = y/—2Q{ri); hence we have 

XB{l-cosT]) (5.14) 



drj 
dl 



with B as above. Equation (5.14) is immediately inte- 
grated, providing 



ctg(77/2) 



XB{z-zq) 



(5.15) 



We can choose zq = 0, so that 77(7r) — 0, and ?/ is anti- 
symmetric with respect to z = 0. In conclusion, the (0,1) 
solution is given by 



?7 = — 2arcctg(Ai3z) 



arccos 



XBz 



%/T+A5 



(5.16) 



this is shown in fig. 4. Note that some care should be 
taken in using appropriate determination of the arcctg 
and arccos functions so that rj is continuous at z = 0. 



C. Comparison with standard Y solitons 

The standard Yakushevich model predictions are re- 
covered for r = A, i.e. for A = 1. The solitons shape is 



evidently very similar to those of standard Yakushevich 
solitons also for A ^^ 1, see figs. 3 and 4. 

In order to compare more precisely the results obtained 
within and without the contact approximation, we recall 
that with the contact approximation the (1,0) and (0,1) 
Yakushevich solitons are given respectively by \12L |20| 



(p{z) ~ 4 arctan [e^^] 
r]{z) — 2 arccos [—Bz/\/Y 



B^ 



,7? = 



(5.17) 



As for the soliton width, which represents the size of 
the transcription bubbles in the Engiander et al. theory, 
this can be determined exactly via (5.9) and (5.15) once 
we decide how this should be measured. That is, the 
width will correspond to z^ — z_ where z± are the points 
at which the angles f or rj differ by their asymptotic value 
(0 or 27r) by less than a given amount A. This is shown 
in fig. 5 for small values of A. 

It should be stressed that the soliton widths vary quite 
slowly with A; thus, the predictions of the model will not 
sensitively depend on the precise value of A. 

More precisely, these are given by the parameter B~^, 
which - see (5.8) - is written as 



B- 



(A/2) ^W{K,r^) ; (5.18) 



again we stress that the curves plotted in fig. 5 have a not 
so steep slope, which shows that the predictions of the 
model do not depend too much on the precise value of 
B. 

In the limit w ^ we have lul = K^6^ and hence 



B- 



b 5 with b 



Kpr^ 



(5.19) 



As clear from fig. 5, dropping the contact approxima- 
tion will cause a widening of the Yakushevich solitons; 
this goes in the right direction since the standard Yaku- 
shevich model produces the right order of magnitude for 
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FIG. 6; Half-width for the (1,0) soliton (upper) and the (0,1) 
soliton (lower) for the Yakushevich model without the contact 
approximation. The curves show the half-width k (in units 
of B~^ = bS) of the solitons as functions of A, with different 
conventions for the measure of the half-width itself, see text: 
A = 7r/20 (continuous curves), A = tt/IO (dotted curves) and 
A = 7r/5 (dashed curves). 



the soliton width but with a too small exact numerical 
value 11]. 

Indeed, it is experimentally known that the transcrip- 
tion bubbles have a width of about 15-20 base pairs; this 
is the size of the region in which the base pairs are open 
and the base sequence can be accessed by the RNA Poly- 
merase. It is not easy to assess what is precisely the angle 
at which base pairs should be considered as open, so we 
have plotted different possibilities in fig. 5. 



VI. CONCLUSIONS AND OUTLOOK 

We have considered the Yakushevich model beyond the 
contact approximation £q -^ 0, focusing on the solitonic 
excitations which - according to Englander et al. ~ are 
supposed to play a functional role in the transcription 
process. 

We have shown that in this case soliton solutions are 
still present, and can be described exactly in analytical 
terms; the simple form obtained within the contact ap- 
proximation is replaced by a more complex expression, 
involving elliptic integrals. However, the qualitative form 
of soliton solutions is little changed, and their width is 
only very moderately increased. 



This shows that the Yakushevich model is actually 
quite robust against changes in £o and dropping of the 
contact approximation; this not really for what concerns 
its mathematical aspects, but rather for what concerns 
its physical features and predictions, in particular in the 
fully nonlinear regime. 

Thus, the first outcome of our work is that one is phys- 
ically quite justified in considering the simplifying con- 
tact approximation in the Yakushevich model, albeit the 
analysis can be performed with the same completeness 
without that approximation. 

Let us mention that other work in progress [5| show 
that by considering a more detailed description - in vari- 
ous ways - of the DNA molecule, one obtains indeed new 
features with respect to simple idealized models as the 
Yakushevich one. 

In this respect, the present work suggests that in ana- 
lyzing these more detailed - and hence more difficult to 
study - models one can in the first instance adopt the 
same kind of approximation adopted by Yakushevich in 
her original study, and focus instead on other features of 
the model. This suggestion will be taken up in forthcom- 
ing work l^]- 



Appendix. Dispersion relations 

In this note we are mainly interested in solitonic ex- 
citations, hence fully nonlinear dynamics. However, the 
study of small amplitude excitations has some interest, 
both per se and in order to emphasize how crucial the 
contact approximation is in this regime. 

Small amplitude dynamics around the equilibrium po- 
sition ■(/)„ = Xn = is described by the linearization of 
(1.10) at (0,0). As Vq, hence V, is non quadratic there, 
these collapse to a pair of identical equations for tp and 
X- In particular, the dispersion relations are now com- 
pletely degenerated, and the intrapair term has no role 
in them. 

This degeneration is removed by introducing in the 
model the "helicoidal" terms mentioned in sect.l 8, lOl 
ri2| . This amounts to introducing in the Lagrangian a 
new potential term 



^^-!^E 



,(1) 



a(2) 



o(2) 
'i+p 



,(1) 



(Al) 
Here p is the half-pitch of the helix in nucleotide units; 
it takes the value p = 5 in B-DNA. 

The introduction of this terms in L entails that a new 
term should be added to the right hand side of the Euler- 
Lagrange equations (1.10); this new term is linear and 
thus is also present in their linearization. 

With standard computations, the new linearized equa- 
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FIG. 7: Dispersion relations for different versions of the Yaku- 
shevich model, (a): The dispersion relations for oj = uj{k) 
as described by (A. 3) for the helicoidal Y model without the 
contact approximation, (b): Dispersion relations for the stan- 
dard (i.e. ^0 = 0) Y model with helicoidal terms, (c): Dis- 
persion relations for the standard Y model without helicoidal 
terms. In all plots, the continuous line refers to the ip branch, 
the dashed one to the x branch. We have used the values 
Ks = 0.13eV/rad^ K^ = 0.009eV/rad^ Q; here k is given 
in units of tt/5, and uj in ps~^. The x branch of the standard 
Y model without helicoidal terms (plot on the right) coincides 
with the degenerate dispersion relations of the Y model with- 
out the contact approximation and without helicoidal terms. 



tions are (note the sign differences in the new terms) 

+ Kh {lpn+p - 2lpn + fpn-p) , 
IXn = Ks iXn+l - 2Xn + Xn-l) + 
- Kh iXn+p + 2X« + Xn-p) ■ 



{A.2) 



8 



Passing to the continuum approximation and Fourier 
transforming via $(a;, t) = /^^ exp[i{kx + ujt)], S(a;, t) = 
Qkuj exp[i{kx + ujt)], the above yield 



OJ. 



{2Ksim- 



u;'= {2Ks/m- 



cos(fc(5)] 
cos(fc(5)] 



mu/i)[i 

i'^Kn/m- 



cos{kp5)\ 
cos{kp5)] 
(A.3) 
for the V' and the x branch respectively (note that here 
the continuous wave number k has the dimension of 
[L]^^, and the lattice spacing 5 sets the space scale; one 
could of course also pass to a dimensionless wave number 
K, — kS). These are the dispersion relations - plotted in 
fig. 5a - for the helicoidal Yakushevich model without the 
contact approximation. Note that a nonzero £o causes 
the presence, even in the helicoidal case, of phonon modes 
(contrary to the ^o = case). 

For comparison purposes, we note that the dispersion 
relations for the standard Yakushevich model are 

cj2 = {2Ks/I)[l - cos{kS)] + {2Kh/I)[l - cos{kpS)] + 

+2Kpr^/I 
uj\ = {2KslI)\\ - cos(fc(5)] -f {2KhlI)\l + co&{kpS)\ 

(A.4) 
(the non-helicoidal case is obtained setting K^ in the 
above); these are plotted in fig. 5b, and in fig. 5c for the 
non-helicoidal (Kh = 0) case. 

Quantities of physical interest are readily evaluated, 
at least numerically, from (A. 3). For Kh ^ 0, the x 
branch has a nontrivial minimum; this is reached for 
k = ko = O.lOlA-i, i.e. for A = Aq = 62A = 18S. It 
is suggestive to remark that if solitons grow out of os- 
cillations triggered by thermal excitations, and assuming 
equipartition of energy between Fourier modes, the larger 
amplitude excitations have a characteristic size which is 
of the order of the size of the transcription bubble, the 
latter being 15 — 20 base pairs [Tj, |2Q||. To the above 
value of fco corresponds a frequency ujq = 0.4ps~^, i.e. 
a period of oscillations Tq ~ 16ps; typical values for Tq 
observed in experiments are in the picoseconds range. 

The speed of phonon excitations Vs corresponds to 
dLu^(k)/dk for fc ^ 0. In the ^o = case, with the values 
of parameters given above, we get Vs — 280m/s, while if 
we drop the contact approximation we get Vg ~ 470m/s. 
These values should be compared with experimental ob- 
servations for the speed of torsional waves, which provide 
values in the range of Vg = 1.6 ± 0.3 Km/s pO| . 
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